In [1]:
library_load <- suppressMessages(
    
    list(
        
        # Seurat 
        library(Seurat), 
        
        # Data 
        library(tidyverse), 
        
        # Reticulate
        library(reticulate), 
        
        # Plotting 
        library(ggplot2)
        
    )
)
In [2]:
random_seed <- 42
set.seed(random_seed)
In [3]:
options(warn=-1)
In [4]:
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
In [5]:
# Source files
source("plotting_global.R")
source("bin/SeuratQC.R")

# SingleRPlot
source("bin/SingleRQC.R")

SCTransform¶

In [6]:
so <- readRDS("data/object/qc.rds")
so$integrate <- paste0(so$treatment, "-", so$sample_rep)
so <- SplitObject(so, split.by="integrate")
so <- so[c("NaCl-Rep2", "NaCl-Rep1", "CpG-Rep2", "CpG-Rep1")]

so <- lapply(so, function(so) {
    
    so <- SCTransform(so, return.only.var.genes=FALSE, variable.features.n=8000, verbose=FALSE)
    so <- RunPCA(so, npcs=100, assay="SCT", verbose=FALSE)
    so <- FindNeighbors(so, dims=1:30, k.param=50, verbose=FALSE)
    so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
    so <- RunUMAP(so, dims=1:100, n.neighbors=50, verbose=FALSE)
    
    return(so)
    
}
                )

saveRDS(so, "data/object/so_sct.rds")
# so <- readRDS("data/object/so_sct.rds")
In [7]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so[[1]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so[[1]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so[[1]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so[[1]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so[[1]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so[[1]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so[[1]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so[[1]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so[[1]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so[[1]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so[[1]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [8]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so[[2]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so[[2]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so[[2]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so[[2]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so[[2]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so[[2]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so[[2]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so[[2]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so[[2]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so[[2]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so[[2]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [9]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so[[3]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so[[3]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so[[3]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so[[3]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so[[3]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so[[3]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so[[3]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so[[3]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so[[3]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so[[3]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so[[3]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [10]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so[[4]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so[[4]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so[[4]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so[[4]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so[[4]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so[[4]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so[[4]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so[[4]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so[[4]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so[[4]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so[[4]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=2000)¶

In [11]:
features_2000 <- SelectIntegrationFeatures(object.list=so, nfeatures=2000)
so <- PrepSCTIntegration(object.list=so, anchor.features=features_2000)
anchors <- FindIntegrationAnchors(object.list=so, normalization.method="SCT", anchor.features=features_2000, verbose=FALSE)
so_int <- IntegrateData(anchorset=anchors, normalization.method="SCT", verbose=FALSE)
In [12]:
so_int <- RunPCA(so_int, npcs=100, assay="integrated", verbose=FALSE)
so_int <- FindNeighbors(so_int, dims=1:30, k.param=50, verbose=FALSE)
so_int <- FindClusters(so_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_int <- RunUMAP(so_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_int, "data/object/so_sct_int_hvg2000.rds")
# so_int <- readRDS("data/object/so_sct_int_hvg2000.rds")
In [13]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=4000)¶

In [14]:
features_4000 <- SelectIntegrationFeatures(object.list=so, nfeatures=4000)
so <- PrepSCTIntegration(object.list=so, anchor.features=features_4000)
anchors <- FindIntegrationAnchors(object.list=so, normalization.method="SCT", anchor.features=features_4000, verbose=FALSE)
so_int <- IntegrateData(anchorset=anchors, normalization.method="SCT", verbose=FALSE)
In [15]:
so_int <- RunPCA(so_int, npcs=100, assay="integrated", verbose=FALSE)
so_int <- FindNeighbors(so_int, dims=1:30, k.param=50, verbose=FALSE)
so_int <- FindClusters(so_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_int <- RunUMAP(so_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_int, "data/object/so_sct_int_hvg4000.rds")
# so_int <- readRDS("data/object/so_sct_int_hvg4000.rds")
In [16]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=6000)¶

In [17]:
features_6000 <- SelectIntegrationFeatures(object.list=so, nfeatures=6000)
so <- PrepSCTIntegration(object.list=so, anchor.features=features_6000)
anchors <- FindIntegrationAnchors(object.list=so, normalization.method="SCT", anchor.features=features_6000, verbose=FALSE)
so_int <- IntegrateData(anchorset=anchors, normalization.method="SCT", verbose=FALSE)
In [18]:
so_int <- RunPCA(so_int, npcs=100, assay="integrated", verbose=FALSE)
so_int <- FindNeighbors(so_int, dims=1:30, k.param=50, verbose=FALSE)
so_int <- FindClusters(so_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_int <- RunUMAP(so_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_int, "data/object/so_sct_int_hvg6000.rds")
# so_int <- readRDS("data/object/so_sct_int_hvg6000.rds")
In [19]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=8000)¶

In [20]:
features_8000 <- SelectIntegrationFeatures(object.list=so, nfeatures=8000)
so <- PrepSCTIntegration(object.list=so, anchor.features=features_8000)
anchors <- FindIntegrationAnchors(object.list=so, normalization.method="SCT", anchor.features=features_8000, verbose=FALSE)
so_int <- IntegrateData(anchorset=anchors, normalization.method="SCT", verbose=FALSE)
In [21]:
so_int <- RunPCA(so_int, npcs=100, assay="integrated", verbose=FALSE)
so_int <- FindNeighbors(so_int, dims=1:30, k.param=50, verbose=FALSE)
so_int <- FindClusters(so_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_int <- RunUMAP(so_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_int, "data/object/so_sct_int_hvg8000.rds")
# so_int <- readRDS("data/object/so_sct_int_hvg8000.rds")
In [22]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Store SCTransform scaled data as h5ad¶

In [23]:
# Store data as h5ad 
adata <- import("anndata", as="adata", convert=FALSE)
pd <- import("pandas", as="pd", convert=FALSE)
np <- import("numpy", as="np", convert=FALSE)

# Re-load SCT object so that it is not PrepSCTIntegration modified
so <- readRDS("data/object/so_sct.rds")
    
# Get SCT scale.data
X_1 <- GetAssayData(so[[1]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_2 <- GetAssayData(so[[2]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_3 <- GetAssayData(so[[3]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_4 <- GetAssayData(so[[4]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()

features <- c(colnames(X_1), colnames(X_2), colnames(X_3), colnames(X_4))
features <- table(features)
features <- features[features==4] %>% names()

X_1 <- X_1[, features]
X_2 <- X_2[, features]
X_3 <- X_3[, features]
X_4 <- X_4[, features]

X <- rbind(X_1, X_2, X_3, X_4)

obs <- rbind(

    so[[1]]@meta.data, 
    so[[2]]@meta.data, 
    so[[3]]@meta.data, 
    so[[4]]@meta.data

)

uns <- list(
 
    hvg_int_2000 = features_2000, 
    hvg_int_4000 = features_4000, 
    hvg_int_6000 = features_6000,
    hvg_int_8000 = features_8000
    
)
    
adata_sct <- adata$AnnData(X=X, obs=obs)
adata_sct$var_names <- colnames(X)
adata_sct$uns <- uns

# Get raw counts
so <- readRDS("data/object/qc.rds")
X_raw <- GetAssayData(so, assay="RNA", slot="counts") %>% as.matrix() %>% t()

adata_raw <- adata$AnnData(X=X_raw, obs=obs)
adata_raw$var_names <- colnames(X_raw)
adata_sct$raw <- adata_raw

adata_sct$write_h5ad("data/object/so_sct.h5ad")

None

SCTransform + regression¶

In [24]:
so <- readRDS("data/object/qc.rds")
so$integrate <- paste0(so$treatment, "-", so$sample_rep)
so <- SplitObject(so, split.by="integrate")
so <- so[c("NaCl-Rep2", "NaCl-Rep1", "CpG-Rep2", "CpG-Rep1")]

so_reg <- lapply(so, function(so) {
    
    so <- SCTransform(so, return.only.var.genes=FALSE, variable.features.n=8000, vars.to.regress=c("msCC_diff_RNA"), verbose=FALSE)
    so <- RunPCA(so, npcs=100, assay="SCT", verbose=FALSE)
    so <- FindNeighbors(so, dims=1:30, k.param=50, verbose=FALSE)
    so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
    so <- RunUMAP(so, dims=1:100, n.neighbors=50, verbose=FALSE)
    
    return(so)
    
}
                )

saveRDS(so_reg, "data/object/so_sct_reg.rds")
# so_reg <- readRDS("data/object/so_sct_reg.rds")
In [25]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg[[1]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg[[1]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg[[1]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg[[1]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg[[1]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg[[1]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg[[1]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg[[1]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg[[1]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg[[1]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg[[1]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [26]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg[[2]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg[[2]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg[[2]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg[[2]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg[[2]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg[[2]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg[[2]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg[[2]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg[[2]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg[[2]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg[[2]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [27]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg[[3]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg[[3]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg[[3]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg[[3]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg[[3]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg[[3]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg[[3]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg[[3]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg[[3]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg[[3]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg[[4]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
In [28]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg[[4]], reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg[[4]], reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg[[4]], reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg[[4]], reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg[[4]], reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg[[4]], reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg[[4]], reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg[[4]], reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg[[4]], reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg[[4]], reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg[[4]], reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=2000)¶

In [29]:
features_2000_reg <- SelectIntegrationFeatures(object.list=so_reg, nfeatures=2000)
so_reg <- PrepSCTIntegration(object.list=so_reg, anchor.features=features_2000_reg)
anchors_reg <- FindIntegrationAnchors(object.list=so_reg, normalization.method="SCT", anchor.features=features_2000_reg, verbose=FALSE)
so_reg_int <- IntegrateData(anchorset=anchors_reg, normalization.method="SCT", verbose=FALSE)
In [30]:
so_reg_int <- RunPCA(so_reg_int, npcs=100, assay="integrated", verbose=FALSE)
so_reg_int <- FindNeighbors(so_reg_int, dims=1:30, k.param=50, verbose=FALSE)
so_reg_int <- FindClusters(so_reg_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_reg_int <- RunUMAP(so_reg_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_reg_int, "data/object/so_sct_reg_int_hvg_2000.rds")
# so_reg_int <- readRDS("data/object/so_sct_reg_int_hvg_2000.rds")
In [31]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=4000)¶

In [32]:
features_4000_reg <- SelectIntegrationFeatures(object.list=so_reg, nfeatures=4000)
so_reg <- PrepSCTIntegration(object.list=so_reg, anchor.features=features_4000_reg)
anchors_reg <- FindIntegrationAnchors(object.list=so_reg, normalization.method="SCT", anchor.features=features_4000_reg, verbose=FALSE)
so_reg_int <- IntegrateData(anchorset=anchors_reg, normalization.method="SCT", verbose=FALSE)
In [33]:
so_reg_int <- RunPCA(so_reg_int, npcs=100, assay="integrated", verbose=FALSE)
so_reg_int <- FindNeighbors(so_reg_int, dims=1:30, k.param=50, verbose=FALSE)
so_reg_int <- FindClusters(so_reg_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_reg_int <- RunUMAP(so_reg_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_reg_int, "data/object/so_sct_reg_int_hvg4000.rds")
# so_reg_int <- readRDS("data/object/so_sct_reg_int_hvg4000.rds")
In [34]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=6000)¶

In [35]:
features_6000_reg <- SelectIntegrationFeatures(object.list=so_reg, nfeatures=6000)
so_reg <- PrepSCTIntegration(object.list=so_reg, anchor.features=features_6000_reg)
anchors_reg <- FindIntegrationAnchors(object.list=so_reg, normalization.method="SCT", anchor.features=features_6000_reg, verbose=FALSE)
so_reg_int <- IntegrateData(anchorset=anchors_reg, normalization.method="SCT", verbose=FALSE)
In [36]:
so_reg_int <- RunPCA(so_reg_int, npcs=100, assay="integrated", verbose=FALSE)
so_reg_int <- FindNeighbors(so_reg_int, dims=1:30, k.param=50, verbose=FALSE)
so_reg_int <- FindClusters(so_reg_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_reg_int <- RunUMAP(so_reg_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_reg_int, "data/object/so_sct_reg_int_hvg6000.rds")
# so_reg_int <- readRDS("data/object/so_sct_reg_int_hvg6000.rds")
In [37]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Integrate (nfeatures=8000)¶

In [38]:
features_8000_reg <- SelectIntegrationFeatures(object.list=so_reg, nfeatures=8000)
so_reg <- PrepSCTIntegration(object.list=so_reg, anchor.features=features_8000_reg)
anchors_reg <- FindIntegrationAnchors(object.list=so_reg, normalization.method="SCT", anchor.features=features_8000_reg, verbose=FALSE)
so_reg_int <- IntegrateData(anchorset=anchors_reg, normalization.method="SCT", verbose=FALSE)
In [39]:
so_reg_int <- RunPCA(so_reg_int, npcs=100, assay="integrated", verbose=FALSE)
so_reg_int <- FindNeighbors(so_reg_int, dims=1:30, k.param=50, verbose=FALSE)
so_reg_int <- FindClusters(so_reg_int, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
so_reg_int <- RunUMAP(so_reg_int, dims=1:100, n.neighbors=50, verbose=FALSE)

saveRDS(so_reg_int, "data/object/so_sct_reg_int_hvg8000.rds")
# so_reg_int <- readRDS("data/object/so_sct_reg_int_hvg8000.rds")
In [40]:
options(repr.plot.width=20, repr.plot.height=15)

dplot_1 <- dplot(so_reg_int, reduction="umap", group_by="seurat_clusters") + ggtitle("Louvain")
dplot_2 <- dplot(so_reg_int, reduction="umap", group_by="treatment") + scale_color_manual(values=color$treatment) + ggtitle("Treatment")
dplot_3 <- dplot(so_reg_int, reduction="umap", group_by="sample_rep") + ggtitle("Replicate")
dplot_4 <- dplot(so_reg_int, reduction="umap", group_by="tissue") + scale_color_manual(values=color$tissue) + ggtitle("Tissue")
dplot_5 <- dplot(so_reg_int, reduction="umap", group_by="cc_phase_class") + scale_color_manual(values=color$cc_phase_class) + ggtitle("Cell cycle")
dplot_6 <- dplot(so_reg_int, reduction="umap", group_by="label_fine_haemosphere") + scale_color_manual(values=color$label_fine_haemosphere) + ggtitle("Haemosphere label")
fplot_1 <- fplot(so_reg_int, reduction="umap", features="pMt_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Mt")
fplot_2 <- fplot(so_reg_int, reduction="umap", features="pHb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Hb")
fplot_3 <- fplot(so_reg_int, reduction="umap", features="pRb_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Percentage Rb")
fplot_4 <- fplot(so_reg_int, reduction="umap", features="nCount_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("UMI")
fplot_5 <- fplot(so_reg_int, reduction="umap", features="nFeature_RNA") + scale_colour_viridis_c(option="inferno") + ggtitle("Feature")

dplot_comb <- dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + dplot_6 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(ncol=4)
dplot_comb

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.

Store SCTransform scaled data as h5ad¶

In [41]:
# Store data as h5ad 
adata <- import("anndata", as="adata", convert=FALSE)
pd <- import("pandas", as="pd", convert=FALSE)
np <- import("numpy", as="np", convert=FALSE)

# Re-load SCT object so that it is not PrepSCTIntegration modified
so <- readRDS("data/object/so_sct_reg.rds")
    
# Get SCT scale.data
X_1 <- GetAssayData(so[[1]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_2 <- GetAssayData(so[[2]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_3 <- GetAssayData(so[[3]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()
X_4 <- GetAssayData(so[[4]], assay="SCT", slot="scale.data") %>% as.matrix() %>% t()

features <- c(colnames(X_1), colnames(X_2), colnames(X_3), colnames(X_4))
features <- table(features)
features <- features[features==4] %>% names()

X_1 <- X_1[, features]
X_2 <- X_2[, features]
X_3 <- X_3[, features]
X_4 <- X_4[, features]

X <- rbind(X_1, X_2, X_3, X_4)

obs <- rbind(

    so[[1]]@meta.data, 
    so[[2]]@meta.data, 
    so[[3]]@meta.data, 
    so[[4]]@meta.data

)

uns <- list(

    hvg_int_2000 = features_2000_reg, 
    hvg_int_4000 = features_4000_reg, 
    hvg_int_6000 = features_6000_reg,
    hvg_int_8000 = features_8000_reg
    
)
    
adata_sct <- adata$AnnData(X=X, obs=obs)
adata_sct$var_names <- colnames(X)
adata_sct$uns <- uns

# Get raw counts
so <- readRDS("data/object/qc.rds")
X_raw <- GetAssayData(so, assay="RNA", slot="counts") %>% as.matrix() %>% t()

adata_raw <- adata$AnnData(X=X_raw, obs=obs)
adata_raw$var_names <- colnames(X_raw)
adata_sct$raw <- adata_raw

adata_sct$write_h5ad("data/object/so_sct_reg.h5ad")

None